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^ Abstract 

In systems with weak ergodicity breaking, the equivalence of time averages and ensemble averages is known to be broken. We study 
here the computation of the power spectrum from realizations of a specific process exhibiting 1 1 noise, the Rebenshtok-Barkai 
model. We show that even the binned power spectrum does not converge in the limit of infinite time, but that instead the resulting 
1—5 value is a random variable stemming from a distribution with finite variance. However, due to the strong correlations in neighboring 
ly-^ frequency bins of the spectrum, the exponent p can be safely estimated by time averages of this type. Analytical calculations are 
,—1 illustrated by numerical simulations. 
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1. Introduction 

In recent years, the study of long range correlations in ex- 
perimental data has moved into the focus of interest. Whereas 
a large amount of empirical results has been obtained using 
detrended fluctuation analysis the same information is, in 
principle, also contained in the power spectrum. Observations 
of nontrivial power law decays of the power spectrum have a 
long tradition and are usually denoted by 1 //-noise, or, more 
precisely, 1 1 with yS =s 1 . Examples include noise in electric 
fesistors, in semiconductor devices (flicker noise), and fluctua- 
tions in geophysical data flSHl- 

\ By weak ergodicity breaking one describes the behavior of 
a system which is ergodic in the sense that a single trajec- 
tory is able to explore completely some invariant component 
of the phase-space (it is indecomposable), but where the invari- 
ant measure is not normalizable, so that the time scales for this 
exploration might diverge (in other words, there is no finite mi- 
croscopic time scale). This concept goes back to Bouchaud Jstl- 
A consequence of this is that time averages computed on a sin- 
gle (infinite) trajectory will not converge to the ensemble mean. 
Instead, the time average is a random variable itself, as a func- 
tion of the initial condition of the trajectory, which is drawn 
from a distribution with non-zero variance. It has been shown 
that Continuous-Time Random Walks |6], but also specific de- 
terministic dynamical systems [7] exhibit this behavior A natu- 
ral consequence of it is the non-existence of a finite correlation 
time, reflected by a power law decay of the auto-correlation 
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function, which translates itself into a 1 //''-behavior of the 
power spectrum. This type of spectrum has been derived for 
a two state process of this type by Margolin and Barkai ^ 
Eq. (33)]. 

As said, the crucial aspect of systems with weak ergodicity 
breaking is that a time average in the limit of an infinite trajec- 
tory does not converge to a sharp value but is a random variable 
with a distribution of finite width. This should also hold for the 
power spectrum. Therefore, it is not evident that an estimation 
of p from a numerically computed power spectrum of a single 
trajectory is reliable. We therefore investigate the estimation of 
the power spectrum from single realizations of a specific pro- 
cess with weak ergodicity breaking. 

A relevant remark is necessary: If one estimates the power 
spectrum by a discrete Fourier transform from a time series with 
sampling interval Af (i.e., uses the periodogram), then the vari- 
ance of the estimator fk for the power contained in the kxh dis- 
crete frequency is f^. However, if the correlations in the time 
domain decay exponentially fast, then the errors in adjacent fre- 
quencies fk and fk±\ are sufficiently weakly correlated such that 
the error of a binned spectrum where one takes averages over 
M adjacent frequencies decays like 1/ Vm [9, Sect. 13.4]. In 
the limit of an infinite time series, the estimation error of the 
binned power spectrum of an ergodic process therefore decays 
to zero, i.e., the spectrum assumes uniquely defined values. In 
the following, we will discuss results for binned power spectra 
of processes with weak ergodicity breaking. The binning im- 
plies that in the analytical calculations we are allowed to ignore 
the time discretization which is present in every numerical data 
analysis. 

The next Sections are devoted to an analytical derivation of 
the fact that the binned power spectrum of a weakly non-ergodic 
process is a random variable whose distribution has a non-zero 
variance. In addition, we compute the correlations between 
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different bins. We will show that because of these correla- 
tions, further binning cannot reduce the uncertainty about the 
true power, and moreover, that the whole uncertainty about the 
power spectrum reduces to a random normalization factor This 
has the important consequence that the power law exponent [3 
can indeed be estimated numerically regardless of the difficul- 
ties to estimate the spectrum itself. 

In Section |5] we illustrate the analytical (asymptotic) results 
by numerical simulations, also opposing time averages to en- 
semble averages. Note that because of the long range correla- 
tions, cutting a long trajectory into pieces and treat these as an 
ensemble is not valid, since this ensemble would not represent 
an independent sample of initial conditions. 

2. The model and basic approach 

Rebenshtok and Barkai introduced the following model for 
a thermodynamic system showing weak ergodicity breaking 
In its simplest form it is characterized by two dis- 
tributions, one waiting time distribution with density 0(f) and 
one distribution of an observable jc e ]R according to the proba- 
bility density k{x). In the model \Qtxo,X\^ ... be i.i.d. random 
variables distributed according to k{x) and let tq,t\,. . . be i.i.d. 
waiting times distributed according to <p{t). The process X{f) 
takes the value xo in the time < f < tq and x\ for the next 
time interval of length ri . In general 



X{t)^Xi forr,_i <r<r,- 
with 



(1) 



(2) 



We assume that the first four moments of k(x) are finite, i.e., 

I^i = J dx x'k(x) < CO for / = 1 , 2, 3, 4. (3) 

Moreover, we assume that the distribution is centered, i.e., jUi = 
0. The waiting time distribution (p(t) should be in the domain of 
normal attraction of an one-sided Levy stable distribution with 
exponent a (0 < a < 1). Therefore, the Laplace transform ^(A) 
of 0(0, 



can be expanded as 

4>iA) = 1 - (v^)" + oiA") as ^ Oh- 



(4) 



(5) 



with the scaling parameter v > 0. 

The Fourier transform of the time series X{t) up to time T is 
given by 



X 



Ft{cS)^ Ate'"'X{t). 



(6) 



The spectrum may be estimated from 
1 



n(T) 



Ft{(jS)Ft(—(jS). 



(7) 



The function n{T) denotes a normalization. Normally, it takes 
the form n{T) - T, but this model requires to take n(T) = T" 
for a non trivial spectrum. For finite time series one would have 
to add correction terms for the fact that the basis for Eq. (|7]| is a 
biased estimate of the correlation function (e.g., see il2l chapter 
8. 1]). These terms vanish for the asymptotic case T ^ oo which 
we are considering here. If we have an ensemble of realizations, 
one can look at the unbinned spectrum 



5ub(w) = Hm 5r(w). 

7"— >oo 



(8) 



Denoting the ensemble average by (■), we are interested in the 
expectation value of the spectrum 



<5ub(w)) = Hm ^{Ft{co)Ft{-uS)) 
and the covariances 

<5ub(Wi)5ub(W2)) 

= lim -4-<Fr(wi)F7-(-wi)Fr(w2)fr(-W2)). 
The limits Eqs. (|9]) and ( fTOb can be rewritten to 



<5ub(w)) = Um — <F,.7-,(w)F,.7-,(-dL»)) 

<5ub(<Wl)5ub('y2)) 

= lim ■\-{F,.Ti(,(Ol)FrT.(,-iLl\)FrT^{iLl2)FrTS-'^2)) 



(9) 



(10) 



(11) 

In the following, it is helpful to define the two- and four point 
correlations 



C2(fl,f2) = <^(fl)^(r2)>, 
C4(fl,f2,f3,f4) = <X(fi)X(f2)X(f3)X(f4)) 

and their (double resp. quadruple) Laplace transforms 

= Jdfi J df2e--^'"-'*^'^C2(fl,f2), 
Ca{Ai,A2, A-i, ^4) = t [C4(fl, f2, f3, f4)] 

d*fe-'"C4(t). 



(12) 



(13) 



(14) 



The double resp. quadruple Laplace transforms of the expres- 
sions in Eq. (fTTI) are 



X[<Fr,(w)Fr,(-w))] 



1 



A\A2 



C^iAx — ia>, A2 + iu)). 



£[{Ft, (wi )Ft,{-co\)Ft, (co2)Ft, (-^2))] 
1 



(15) 



A1A2A3A4 



C4(Ai — iciJi , A2 + ia>i , /I3 - 10)2, A4 + ia>2)- 



The limits in Eq. (fTTl i can also be expressed in Laplace space 
using a formulation following the multidimensional Tauberian 
theorem by Drozhzhinov and Zav'jalov lll3Lll4tl as 



Lf^co "J 

- Um — — C2(^/li - !w,^/l2 + 

X [ lim — <F,.7-, ( W 1 )Fr7-, ( - W 1 )f ^r, (c^2)FrTi (-m))] 



■2a 



lim 

f^0+ /ii/i2/l3/l4 



C4( ^Ai - icoi, (A2 + ia>i , 
- ico2, (A4 + 10)2) ■ 



(16) 



(17) 



In the next section, we will show how one can calculate the 
Laplace transforms C„{A) of the multi point correlations Cnit). 

However, when only one time series is given, taking 5r(w) 
from Eq. O as estimate of the spectrum will lead to unsatisfac- 
tory results as the value of S ricL)) will fluctuate except for some 
pathological cases. One often relies to the concept of binning. 
Instead of considering a single frequency co, one averages all 
available frequencies in an given interval [oj - ^Aw, a> + \Aol>]. 
For a time series of length T the Fourier transform returns val- 
ues at frequencies which are evenly spaced with distance y- 
For large T the sum can be approximated by an integral and we 
can consider as observable for the binned spectrum 



St(oj, Aw) 



1 p 



do)' St (co')- 



The spectrum can be estimated from 
5 bin (t^, Aw) = Hm57(w, Aw). 



(18) 



(19) 



Analytically, it is helpful to let the bin size go to zero as a last 
step 



Sbmi(^)- lim 5(w, Aw). 



(20) 



It is important to take the last limit at the very end. In most of 
the important cases, the value of Suni<^) will converge almost 
surely to the value of the spectrum. 

In this paper, we want to consider the situation where the 
stochastic process X(t) is the model by Rebenshtok and Barkai. 
The estimation of the spectrum is connected with averaging 
over the time series, therefore the question arises how the weak 
ergodicity breaking affects the observable 5bin(w). Does the 
value of 5bin(w) converge to a single value or does it converge 
to a non trivial probability distribution? It will turn out that the 
latter is the case. Unlike the observable discussed by Reben- 
shtok and Barkai 1 10, 11], it seems not to be possible to deter- 
mine directly the probability distribution of Sun((^)- But one 
can obtain important information from the moments (5'bin(w)) 
and <5bin(wi)5bin(w2)). 

The first step is to look at {S (w, Aw)) with 2|w| > Aw: 



(21) 



We proceed similarly to the unbinned case and formulate the 
limit in Laplace space (note that the Laplace transform and 
the w-integral commute by Fubini's theorem and the estimate 
(w)Fj--,(-w))| < fJ.2T\T2 where fi2 is defined in Eq. (O) 



1 r 

r^o^ Aw J„ 

Aw 



dw' {FrT,{(Jj')F,-T2i-^')) 



Aa> 



(22) 



- lim 



dw' 



C2(^Ai-ioj',^A2 + ioj'). 



Similarly, (5(wi, Awi)5(w2, AW2)) with 2|wi| > Awi and 
2|w2| > Aw2 can be obtained from 



<5bin(<Wl , Awi)5bm(W2, AW2)) 



r^oo Awi 



1 r^OJi + -^/\OJi p< 

- — ' dw; 

|AW2 Jaj,-Uoj, J 01 



■aj2 + h AtL»2 



dwo 



{F,.T,{0j'l)FrT2{-<^'l)FrTi{(^2)FrT4i-'^'2)) 



(23) 



r3=7'4=i 



One gets 



lim I dw 

c^oo AwiAw2 J„,-i A«, 



'0)1 +^ Att>| 



= lim 



'1 , dw^ 

kJ OJ2 - \ Ati)2 

{F,T, (w'l )F,r, (- w'l )F,T, iaj'2)Frn(-ca2)} 
I dw'i I dw2 

Jw,-iAtJi Jmz-^Am, AiA2y 



^^0+AwiAw2 Jaj,-iAaJt Jc02-^_Ac02 A1A2A3A4 

C^i^Ai - , ^A2 + ioj'i , ^At, — /wj, f /i4 + iWj). 



(24) 



Therefore the problem can be reduced to the determination of 
the Laplace transforms C„{A) of the n point correlations func- 
tions. We describe in the next section how to obtain them. 



3. A diagrammatic approach to the multi point correlation 
functions 



In III5II a diagrammatic method was introduced how to de- 
termine the joint probability distributions and multi point cor- 
relation functions of a continuous-time random walk. In this 
section we are adapting this method to determine C„(A). For 
any set of natural numbers q\, . . .,q„ € No, one defines 



(25) 



with the indicator function (A being any subset A c M) 



(5 bin (t'^, Aw)) = Hm -— 
Aw 



■w+^Aoj 



dw' <Fr(w')Fr(-w')). 



a)-4 AtJ 



1a(0 = 



|1 forfeA, 
[0 forf^A. 



(26) 



The C„[q](t) can be considered as partial correlations which 
only contribute, when the time f,- is in the {qi + l)th waiting time 
(for / = !,...,«). Therefore, 



C„(t) = <X(fi)---X(f„)> 

CO 

Define the contribution of one step as 

i:qi=j i:qi>i 

One obtains by induction 

,0)r 



(27) 



•■■ii<j 



C„[q](t)= ★;7;,^'[q](t) 



(28) 



(29) 



where ★ denotes the convolution in the time variable f i , . . . , f„ 
and /■ runs from to any natural number larger than 
max(^i,...,^„). 

Transforming into Laplace space gives 



(30) 



Using the following notations from 01511 (suppressing the de- 
pendence on q) 

V; = {/ : = /), Ej^{i:qi<j}, hj^{i:qi> j). (31) 

The set V^- ("vertex") contains the indices which are at position 
j, the set Ej ("earlier") contains the indices which are before 
that position and the set Lj ("later") contains the indices which 
are after that position. With the notation 



A, 



(32) 



one gets 



^ ieVy ieEi ieh: 



2 (-l)l'l0(AyuL,). 



(33) 



Here, the power set of V; is denoted by PCVj) and |7| is the 
number of elements of J. Therefore, the last sum goes over all 
subsets J of indices of Yj. 

As in f]3] we introduce diagrams which denote the relative 
ordering of the ^,s. A generic diagram is given in Fig.lT] The 
main idea behind this consists in grouping all steps / with the 
same V,-, E/ and L,-. If V/ = {) then there can be several suc- 
cessive steps with the same E,- and L,-. These successive steps 
are combined to a horizontal Une. If V, ^ {) then there can 
be only one step with these sets of indices. This is represented 




Lo 



Figure 1 : A generic diagram 



a) 



2 , 

r 



b) 4 3 2 



4 3 




Figure 2: Tlie types of diagrams relevant for a) the two point congelation and b) 
four point correlation 



by a vertex with the indices V, leaving. Therefore the diagram 
Fig-fflrepresents all q with the following properties: First there 
is any number of steps / with E,- = {), V/ - {) and L/ - Lq. 
A single step with E/ - {), V,- - Vi and L,- = Li follows. At 
next is any number of steps / with E; = Vi, V,- = {) and L,- = Li. 
The rest in interpreted analogously. For each n point correlation 
there is only a finite number of diagrams and summing over all 
q represented by a diagram turns out to be easy. It turns out 
that one can associate to each part of the diagram a factor (for 
a more comprehensive derivation of these results, we refer to 
(l5]). One obtains the following rules (where yu, is defined in 
Eq. ©) 



^vertex /' 



n, 



^ (-l)l'l0(AyuL,), 



riir 



CO ^ 



(34) 



mo 



Since we assumed that yUi = 0, we do not need to consider 
diagrams which contain vertices with only one leaving line. 
Therefore a single diagram remains for C2(M, ^^2) which drawn 
in Fig. |2^. This gives 



C2(/ll,/t2) - yiineX=(l,2)yvertex,V={l,2),L={) 

1^2 1 - 0(/ll) - 4>(^2) + 4>(A + ^2) 



(35) 



^1^2 1 - cj)(Ai + ^2) 

For the four point correlation two types of diagrams are rele- 



4 



vant, which are drawn in Fig.|2j3. The contributions are 

DliA) - 'yime,L=|I,2,3,4)Tvertex,V=|3,4),L=(l,2) 
X Tline,L=|I,2)Tvertex,V=(l,2|,L=() 
_ f^l ^(A(i,2|) - 0(A|i,2,3)) - 0(A|1,2,4|) + ^(A|i,2,3,4)) 
/Il/l2/l3^4 1 - ^(A|i_2,3,4)) 

I - 4>(Ai) - 0(A2) + ^(Ai + A2) 



1 



+ ^2) 



(36) 



The case Q 9^ appears if and only if the waiting times t, can 
take only the values 0,^,^, Therefore, more generally, 

^(z) = 0(z + icj) for Re(z) > and w e D.Z. (44) 



In this paper we concentrate on the spectrum for frequencies 
ll) i 0.11. The frequencies oj € QZ behave differently, e.g., 
the behavior of 5ub(0) is completely described by the result of 
Rebenshtok and Barkai ifiolfllll . Now, we can proceed to cal- 
culate the Umit Eq. ( fTTT i by looking at the terms Eqs 
(|37] ) with the parameters given in Eq. ( l38T l. We have 



Dn(A) -7line,L=ll,2,3,4|Tvenex,y=|l,2,3,4|,L={l 

fl4 EjeP((l,2.3.4|)(-l)''^'0(Aj) 
A1A2A2A4 1 - ^(Aji_2,3,4|) 



(37) 



The function C4(A) is obtained by summing over all possible 
combinations of indices at the vertices, i.e.. 



(38) 



C4(i) ^Di(Ai, A2, A3, A4) + DiiAi,A3, A2, A4) 

+ Di(Au A4, A2, Ai) + D^{A2, A-i,Ay,A4) 
+ Uii^i, M, + DiiAi, A4, Au A2) 

+ Dii(Ai,A2,A3,A4). 

4. Evaluating the spectral observables 

4.1. The unbinned observables 

We can use the results of the last section to determine the 
spectral observables. Using the expansion Eq. (|5]l and the two 
point correlations Eq. ( l35T l to calculate the limit Eq. ( fT6b 



lim (^Dii^Ai - icij] , (A2 + ia)\ , (A3 - iu>2, (A4 + ia>2) 



+ Dii^A^ - ia>2, (A4 + ia>2, (Ai — ico] , (A2 + iu){)^ 



A" 

(1,2,3,41 



1 1 
+ 



A"" 

v'Ni,2) 



A"" 

13,4) y 



<5ub(wi))<5ub(W2)). 



(45) 



Using the fact 0(/dL»i)lnz(wi -012) - 4)(icj2)lnz(<^\ -<^2), yields 
similarly 

lim (Dii^Ai - icL>i,(A4 + iu2, (A2 + ia)\,(A3 - ico2) 

+ Di{(A2 + , (A3 - i(i)2, (Ai - icL>i , (A4 + 10^2)) 



= r(i +ay 



A" 

(1,2,3,41 



1 1 
+ 



A" 

V 12,3) 



A"* 



<5ub(Wi))<5ub(W2)) 



X lnz(wi - W2) 



(46) 



and 



Um C2(<r^i - iLo,CA2 + iu) 

f^o+ A1A2 



1 



//2 2 - (j){iu) - 4>{-iu) 



AiA2{Ai + A2Y V" aP- 
The Laplace transform can be inverted by using 



min(ri,r2)'' 



1 



r(l + a) 
Therefore, one has 

1^2 



A^A2{Ai + A2Y 



<5ub(w)) = 



2 - 4>{i(j->) - 0(-;w) 



v«r(l -Ha) 
with the behavior near u 



0+ 



<5ub(w)) = 



2^12 cos (fa) 



1 



v'T(l +a) \cjj\ 



l-a 



+ O 



,\l-a 



(39) 



(40) 



(41) 



(42) 



This last equation shows that the spectrum has a 1 // " form 
near the origin. 

Following lfl6[ Lemma XV. 1.3], there exists an f2 > such 
that 



^(/<jj) — 1 



w e QZ. 



(43) 



lim {p\(C,A\ - /wi, f/}3 - /w2, (A2 H- i(L>i,(A4 + 10)2) 



-I- Di{(A2 + icJi , (A4 + icj2, (Ai - icL)\ , (A3 — /W2)) 



r(l +aY 



A" A" 

^^1,2,3,41 V |1,3) 



1 1 

+ 



A"" 

|2,4)y 



<5ub(Wi))<5ub(W2)) 



X lnz(wi -I- (JJ2). 



Finally 



(47) 



Hm C"Dii((Ai - ioJi , (A4 + ia)2, (A2 + ico^ , (A3 - ia)2) = 0. (48) 



In Appendix A we derive the following Laplace transform 



min(ri,r2, T3, T4f m:m{TuT2f 

^ , min(ri,r2,r3,r4) ^ 

X F{a, -a; \ + a, ; ) 



min(ri,r2) 



(49) 



1 



^1/12^3^4 Aj\_2,3,4| ^(1,2) 



with the hypergeometric function F{a, -a\ I + a\ x). As 



F(a, -a; 1 H- or; 1) = 



r(l + 2a) 



(50) 



we obtain for the limit Eq. (fTTl i 



and accordingly 



r(l +af 

<5ub(Wl)5ub(W2)) =2— r— -<5ub(wi))<5ub(W2)) 

r( 1 + 2a) 

X (1 + lnz(wi - 0)2) + lnz(wi + wi)) . 

(51) 

Eq. dSTI ) will play an important role in determining the prop- 
erties of the binned case. In particular, the last equation deter- 
mines the variance of S^hU^) (for lu i OZ) 



<5bm(w)> = <5ub(w)>. 



(59) 



Var[5ub(w)] = yub(a)<5ub(w))2 
with the function 

r(i + af 

Vub(«)^4 

r(i -I- /ff) 



(52) 



(53) 



As already mentioned in the introduction, even in the ergodic 
limit a ^ 1, the value yub(l) = 1 is larger than zero and there- 
fore the observable 5 ub('^) does fluctuate. 

The correlation coeflicient between 5ub(wi) and5ub(w2) (for 
2(jL>i t QZ and 2(1)2 i i^Z; here Q denotes the periodicity of 
^{X), see Eq. (|43]l) 



p[5ub(Wi),5ub(W2)] = 



1 1 for wi - W2 e 
\R{a) for tL)\- (jj-ii Q.Z 



with 



2r(l H- a)2 - r(l + 2a) 



4r(l H-a)2 -r(l +2a)' 
In many cases we have Q = 0, i.e., 

p[5ub(<wi),5ub(w2)] = R(a) for 1^1 1 1^2! . 



(54) 



(55) 



(56) 



In the ergodic limit, we have R(l) - 0. Therefore, the obser- 
vations 5ub(<yi) and Subi^^i) of the spectrum at two difi'erent 
frequencies cji and (l)2 are uncorrelated. 

4.2. The binned observables 

Here, we are going to use the results from the unbinned case 
to evaluate the limits Eqs. (l22l l and (l24l i. The limits on the right 
hand sides of these equations can be commuted with the fre- 
quency integrals describing the binning. As this statement is 
not obvious, we are giving an argument. 

First, we are considering (5'bin(<J^, Aw)) with 2\u)\ > Aw. 
Then we have the simple estimate 



1 



1 - mM + ^2)) 



(57) 



The commutativity follows by dominated convergence. There- 
fore, we get 



<5bin(w, Aw)) 



dw' <5ub(w')) 



(58) 



It is not surprising that the binned estimate yields also the cor- 
rect expectation value for an ensemble of realizations. 

The argument for {S bin(<wi , Awi )5bin(w2, AW2)) (with 2|wi | > 
Awi and 2|w2| > AW2) needs more work. First notice that we 
have the following inequality 



\1-4>{A + iLo)\ > 1 - \(f>(A + /w)| > 1 - 4>{A). 



(60) 



This gives the following estimates (replace the tilded variables 
by the corresponding permutations of the /l,s and w,s) 

+ ia>\ , ^^2 + i<^2, + ^^4 + /w4)| 
^1 4 4 (61) 

" IW1W2W3W4I 1 - ^(^A(i,2,3,4l) 1 - ^(^A(i,2)) 

and 

|Du(^/Ii — ia>\ , ^/l2 + iu)\ , I^At, — ia>2, (A^ + ia>2)\ 

^ 1 16 (62) 

~ W^w^ 1 -0(^A(i,2,3,4|) 

By dominated convergence and Eq. ( BTT l. this leads to 



{Sbm(a>l , Awi)5bm(W2, AW2)) 

— — dw'i dw2<5ub(w'i)5ub(w2)) 

AW1AW2 Ja„-lAoJt Joj2-^_Aa,2 

'^^^'^"^ I^^J^dw'i J"dW2<5bm(^))<5bin(W2)) 



= 2- 



r(l -I- 2a) Awi 



(63) 



or shorter 



r(l -H a)2 

{Sbm(LL>l)Shin(CL>2)) = 2— —{Shm(CL>l)){Shin(m)) (64) 

r( 1 -H 2a) 

for all wi , W2 0. The additional factor present in Eq. ( fSTT l for 
the values wi - W2 e Q.Z and wi -i- W2 e Q.Z has vanished due 
to the binning. 

The variance of Sbm((^) turns out to be 



Var[5b,n(w)] = Vu„ia){S unioj)}' 
with the function 
V . ^ .r(l+a)2 



(65) 



(66) 



For < ff < 1, the value of Vbin(Q') is positive, i.e., the vari- 
ance of Sbm((^) does not vanish and the observable 5'bin(<J^) is a 
proper probability distribution. ybin(Q;) describes the part of the 
fluctuation of the spectrum which is connected with the weak 
ergodicity breaking, i.e., the part that remains after binning. 
Correspondingly, for the ergodic limit a — > 1, the function 
Vbinia) vanishes (y(l) = 0) and the observable 5bin(w) does 
not fluctuate. 



In general, the whole spectrum is estimated. The question 
arises, how the weak ergodicity breaking affect different fre- 
quencies. Is it possible that one would estimate a wrong ex- 
ponent y6 of the 1//'^ noise? This can be answered by using 
Eq. ( |64] | to calculate the correlation coefficient for all < a < 1 
and (D[,aj2 + 0: 

p[5bm(wi),5bin(W2)] = 1. (67) 

This implies that 5bin(<J^i) and 5bin(w2) are Unearly coupled 



(5bin(Wl)) 



(68) 



Therefore, as soon as one value of the spectrum is evaluated, 
the other values follows from this one. In other words, the weak 
ergodicity breaking does only affect a random prefactor which 
is the same for all frequencies. 

5. Numerical Verification 

In order to verify the presented analytical results we have per- 
formed numerical realizations of the process described above. 
To this purpose we have generated time series according to the 
following procedure: 

• Pick random numbers, r,, from a totally asymmetric Levy 
a-stable distribution with a given a. These numbers rep- 
resent the waiting times. The random numbers were gen- 
erated by means of the procedure introduced in iITtIi with 
the scaling parameter set to one. 

• Pick additional random numbers, Xh from a normal distri- 
bution. These numbers are the values of the process during 
the respective time intervals. 

• The values Xi are written into the entries of an array. Ay, 
of total length L, for ; e [Af 2^ A?X;j.^i Ti[. In all 
simulations presented here we chose Af = 0.25, defining 
the minimal resolution of the process. 

The unbinned spectra, 5iib(w), can than be determined straight 
forwardly by means of the fast Fourier transform of the array 
A. Of course, the spectra obtained for these processes are dis- 
crete. In order to obtain the binned spectra, 5'bin('^, Aw), we av- 
erage 5ub over various consecutive frequencies and assign the 
obtained value to the average of these frequencies. 

In Fig.|3]we plot the unbinned spectrum of a typical realiza- 
tion of a process of length L = 2^' for a - 0.5 together with 
the binned spectra for two choices of the binning. Aw = 2"'^7r, 
and Aw = 2"'"7r, respectively. It should be observed how the 
spectra get smoother when increasing the bin size, while the 
average slope of the spectra is not changed by the binning. De- 
viations from the power law behavior at small and large fre- 
quencies stem from the finite size of array A and the minimal 
resolution of the process, respectively. 

The inset of Fig. |3] shows the binned spectrum with Aw = 
for various realizations of the processes with length L - 
2^' for a - 0.5. Note that, while the exponent characterizing 




Figure 3: Power spectra corresponding to a single realization of the waiting 
time process are plotted for a = 0.5 for the unbinned case as well as for two 
different choices of the bin size. In the inset we show spectra corresponding to 
different realizations of the same process for a fixed choice of the binning for 
a = 0.5. Note that the curves exhibit a shift with respect to one another. 
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Figure 4: The averaged unbinned spectra corresponding to various choices of a 
are plotted. Note that the spectra show an algebraic decay over various decades 
oio) with an exponent /? depending on a. The curves have been shifted for better 
visibiUty. In the inset we compare the theoretical prediction for the exponent /) 
(solid line) to the numerical findings (dots). 



the decay of the spectrum does not depend on the realization, 
the normaUzation factor of the spectrum does. We would like 
to emphasize that the normalization factor will neither converge 
by increasing the length of the process nor by choosing larger 
bins. This probabilistic property of the normalization factor is 
a consequence of the weak ergodicity breaking. 

To verify our predictions with respect to the dependence of 
the spectrum on a, we have considered the unbinned spectrum 
averaged over 10000 realizations of the process, {5ub(w)). In 
Fig. |4]we plot (5ub(w)) for processes of length L = 2^' and 
a - 0.1,0.4,0.7,0.9 (bottom to top), respectively. The curves 
have been shifted with respect to one another for better visibil- 
ity. Note that the spectra show a ~ 1 / w'' behavior over various 
decades with an an exponent yS that decreases with increasing a. 
In the inset of Fig. |4]we compare our prediction for the value of 
the exponent /3{a) - 2 - a (solid line), see Eq. (l42T i. to the val- 
ues obtained numerically by means of measuring the slope of 
the spectra in log-log scale at the intermediate region. A fairly 
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Figure 5: The theoretical predictions for V(a) are plotted for the unbinned (.solid 
line) and binned (dashed line) case. Numerical results for the unbinned and 
the binned case are represented by the square and circular dots, respectively. 
In the binned case, various choices of process length have been considered, 
L = 2'^, 2'^, 2", 2^^ (top to bottom). The inset shows the proportionality factor 
between the spectrum and the square root of the variance for the unbinned case, 
^ub(<^) = Var[5ub('<J)]/(5ub(<^))^. for <^ = 0.5. Note that it is constant over 
various decades. 



good agreement is observed. 

Let us now check the predictions with respect to the variance 
of the spectra. To calculate the variance of the unbinned spectra 
we have considered processes of length L - 2^' and averaging 
was carried out over 10000 realizations of the process. In accor- 
dance with Eq. (l52b it is found that the variance is proportional 
to the squared expectation value of the spectrum over various 
decades with deviations at very high and very small frequen- 
cies only (see inset of Fig.|5]). To determine the proportionality 
factor we have averaged this factor over the intermediate fre- 
quency region. In Fig.|5]the theoretical prediction Eq. ( |53] ) for 
the proportionality factor Vuh(a) (solid line) is compared to the 
numerical results (square dots), for various choices of a. We 
find a reasonable correspondence with deviations at very small 
a. However, these deviations become smaller when considering 
longer process lengths. 

To verify the predictions for the asymptotic behavior of 
binned spectra we have considered processes of various lengths, 
L - 2'-', 2'^, 2'^, 2^-', while the size of the bins in frequency 
space was kept fixed at Aw = l^'^n. Note that, by augmenting 
the process length while keeping the bin size fixed, the number 
of frequencies of the discrete spectrum falling into a specific 
bin increases. Again, the averaging was carried out over 10000 
realizations of each process. In Fig. |5]the theoretical predic- 
tion for ybin(Q') for the binned case (dashed line), see Eq. ( |66] |. 
is compared to the numerical data obtained for the binned pro- 
cesses (circular dots). It should be observed that, as the length 
of the process is increased (top to bottom), the values of Vbin de- 
part from the predictions of the unbinned case (solid line) and 
approach the predictions for the binned case (dashed line). 

Next, we have verified our predictions with respect to the 
correlation coefficient of the spectra. For the unbinned case we 
have again considered processes of length L = 2^' and aver- 
aged over 10000 realizations. In Fig. |6l the correlation coef- 



Figure 6: plS uhiaii), S ^^(012)] is plotted for a>i = 2n and various choices of a. 
The curves shift upwards as a is decreased. The peak appears at aij = li>2- For 
very small a the peak extends over the whole frequency domain. 
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Figure 7: The theoretical predictions for P|tj| |* | (a) are plotted for the un- 
binned (sohd line) and binned (dashed line) case. Numerical results for the 
unbinned and the binned case are represented by the square and circular' dots, 
respectively. In the binned case, various choices of process length have been 
considered, L = 2'^ 2'^, 2", 2",223 (bottom to top). 



ficient of the unbinned spectrum, p[5ub(wi), 5ub('JJ2)], is plot- 
ted for the fixed frequency (jJ\ - 2n and various choices of a. 
Note that the correlation coefficient for not too small a is in- 
deed flat (apart from the trivial peak at (±i\ - 0)2) and that its 
value increases for decreasing a. For very small values of a the 
peak around <jJ\ - u)2 gets extended over the whole frequency 
domain but this eff'ect can be remedied by considering longer 
process lengths. 

In Fig. Q we compare the predicted value for the correlation 
coefficient in the unbinned case (solid line), p|oj,|#Im,|(Q') = R{a), 
see Eq. dSSl l. with the numerically obtained data (squared dots). 
For the numerical data the correlation coefficient was averaged 
over 2'* choices of the pairs wi and W2. Reasonable agree- 
ment between theory and numerical data is found with devia- 
tions at small a values stemming from the broadening of the 
peak around \ - |w2l- 

Finally, we repeated the analysis of the correlation coeffi- 
cient for the binned spectra. Again, we considered processes 
of various lengths, L - 2'-', 2'^, 2'^, 2'^, 2^^, while keeping 



the bin size fixed in frequency space with Aw 



In 



Fig- 13 the theoretical prediction for the correlation coefficient, 
P|w,|^|w2|(a') - 1, see Eq. (l67l l. is plotted as dotted line while the 
numerical data is represented by the circular dots. As the length 
of the process is increased, the numerical values for p\a,,\^\a>^\ 
depart from the theoretical prediction for the unbinned case 
(solid line) and approach the predicted line for the binned case 
(dashed line). 



6. Summary 

In this paper, we have determined the spectral properties of a 
model by Rebenshtok and Barkai which shows weak ergodicity 
breaking. The analytical results were verified by numerical sim- 
ulations. Near the origin, the spectrum shows a typical l/f^^" 
behavior with a characterizes the tail behavior of the waiting 
time. Using a single time series will result in fluctuations of 
the unbinned spectrum which are also seen in the ergodic case. 
Therefore, one commonly uses binning to determine reliable 
values. While the fluctuations of the binned spectral observ- 
able vanish in the ergodic case, this is not the case for weak 
ergodicity breaking. However, the fluctuation does only aff'ect 
a common prefactor for the whole spectrum. Consequently, the 
measurement of the exponent /? of the l/f^ behavior is not hin- 
dered by the weak ergodicity breaking. 

Recently, the emergence of universal fluctuation for pro- 
cesses with weak ergodicity breaking has been discussed by 
He et al. |18], Sokolov et al. |19] and Esposito et al. li20J. It 
has been shown that for several time averaged observables as 
the average mean square displacement the fluctuations are de- 
scribed by a Mittag-Leffler distribution. As can be infered from 
Eqs. (|65] | and (|66] |, the first two moments of the binned spec- 
trum are also in agreement with the assumption of a Mittag- 
Leffler distribution. However, it remains to be seen in a future 
work if the binned spectrum does also belong to this universal- 
ity class. 

The analytical results for the binned spectrum are valid in the 
hmit of infinite time series length. In a numerical reahzation 
of this, binning means averaging only over a finite number of 
discrete frequencies due to the finite process time. It is therefore 
a nontrivial observation that these binned spectra converge for 
numerical manageable process lengths towards the theoretical 
prediction. 
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Appendix A. Calculation of the quadruple Laplace trans- 
form 

In this appendix, we derive the Laplace transform Eq. ( |49] l. 
We want to mention that it is possible to obtain and motivate this 



result by the more general methods discussed in the derivation 
of Eq. (71) in ifisll . First consider the function 



g(t) = 6{n - tixti - hr-'6{t2 - h)6{h - f4)fr' 



(A.l) 



with the Dirac delta d{t) and the Heaviside step function 9{t). 
Its quadruple Laplace transform can be directly calculated: 



T{\+ay 1 1 

" ^^(1,2,3,41 {1,2| 



(A.2) 



Integration yields 



dfi dt2 dti df4g(t) 

Jo Jo Jo Jo 

1 /^min(ri,7'2,r3,7'4) 

= - df (min(ri,r2)-f)"f''"^ 

a Jo 

= J- min(r, , 72)" min(ri , Tz, Ts, 74)" 



(A.3) 



X F ff, -a; 1 -I- ff; 



min(ri,r2, Ts, r4) 
min(ri,r2) 

Eq. ( |49] l follows from the properties of the Laplace transform 
with respect to integration. 
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